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We describe generalized time-dependent 
mean-field equations for partially 
condensed samples of trapped and eva- 
poratively cooled atoms. These equations 
give a way of investigating the various order 
parameters that may be present as well 
as the existence of a mean value of the 
field due to condensed atoms. Our 
approach provides us with a closed system 
of self-consistent equations for the order 
parameters present. The equations we 
derive are shown to reduce to other 
treatments in the literature in various limits. 
We also show how the equation of 
motion method allows us to construct a 
formalism that can handle the evolution 
of these mean fields due to two-loop 
kinetics. 
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1. Introduction to Generalized 
Mean Fields 

A great deal of interest has been generated in Bose- 
Einstein Condensation (BEC) by the recent experimen- 
tal observations of this phenomenon in evaporatively 
cooled alkah gases [1], [2], [3]. These experiments have, 
in particular, given great impetus to the study of in- 
homogeneous Bose gases. There are, of course, a fair 
number of treatments of the Bose gas in the literature, 
which make use of the various tools of the many body 
trade and a thorough account of the special features of 
the physics of homogeneous Bose-condensed liquids is 



given in [4]. Most treatments deal with Bose-condensed 
systems that are static and either in, or close to, thermal 
equilibrium. The majority of studies have also, until 
recently, been restricted to homogeneous gases, with the 
notable exceptions of the work of Pitaevskii [5] and 
Fetter [6] . However, due to the fact that BEC has been 
experimentally achieved in a trap, there are several more 
recent studies that have addressed the inhomogeneous 
case. In this paper we shall present an equation of 
motion approach [7] that is well-suited to handling both 
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time-dependent and non-equilibrium cases for a trapped 
(i.e., inhomogeneous) gas. We shall also show how this 
approach reduces to the conventional T ^ theory in 
the appropriate limits. 

At non-zero temperatures the atomic assembly does, 
of course, consist of excitations above the condensate 
which have to be taken properly into account. Perhaps 
the best known way to do this is based on the finite- 
temperature Hartree-Fock-Bogoliubov (HFB) equations 
[8]-[ll]. These equations include the effects of the 
mean value of the atomic field, the effects of thermally 
excited atoms, as well as the so-called anomalous 
averages of atomic operators [12]. These anomalous 
averages represent strong pair correlations between the 
particles, that can arise when the effective interactions 
are attractive. The HFB equations can be derived in a 
variety of ways, e.g., via Green functions, variational 
principles [8], [11] and through an equation of motion 
approach. An important discussion of the HFB theory 
has been recently given by Griffin [13], who describes 
how the behavior of an inhomogeneous Bose-condensed 
gas can be obtained for the full temperature range within 
the context of the Popov approximation [14] (explained 
below). 

In this paper we present what could be considered a 
generalization of his approach to include the possibility 
of yet further order parameters of the type discussed by 
Laloe [15]. In doing this, we also lay the ground for a 
fully non-equilibrium technique. We have chosen to 
work with the equation of motion method, as opposed 
to other approaches like the non-equilibrium Green 
function technique (Keldysh Method) [16], [17], as it is 
somewhat easier to produce kinetic equations and also 
to give a link to simpler equilibrium theories. Further- 
more, the equation of motion approach makes it easier to 
look at "one-loop" effects away from equilibrium with 
arbitrary initial conditions. We emphasize that our 
techniques are aimed at the dilute gases that are being 
produced by the evaporative cooling method [18]. For 
those, a kinetic equation approach should be a very 
accurate and convenient technique. It is also possible to 
set up efficient simulation techniques based on Monte- 
Carlo wavefunction methods for a generalized density 
matrix approach. 

We still need, however, to construct a proper self- 
consistent starting point for the kinetics and employ a 
consistent decoupling approximation. We believe this to 
be equivalent to the approach advocated in [19] for the 
calculation of the time variation of the superconducting 
order parameter in a Fermi liquid. Central to such an 
approach lies the idea that the main effects of the 
excitations can be described through the mean fields 
("one-loop" diagrams), whereas all higher order 
correlation effects can be summed up into kinetic terms 



terms (i.e., "two-loop" effects). We shall not be able to 
give a full account of the kinetic terms here, but we will 
show clearly how they enter our formulation and how 
they can be computed. The full details of the kinetic 
terms will be given elsewhere [20]. 

We start our analysis by considering a group of atoms 
trapped by an external field. The spatial and temporal 
behavior can both be fully described, in second quan- 
tization, in terms of the boson field operator ^(r, t). 
We will first outline our approach in terms of ^(r, t), 
before giving our results for the occupation number 
representation in the subsequent sections, the latter 
treatment being better suited to a computational 
approach to the problem. For the dilute gases under 
consideration, we assume that the dominant interatomic 
interactions occur pairwise. The hamiltonian of the 
system can be thus written in the form 



H= l^PXr, 



t) 



2m 



-^Vtmp (r,t) 



^(r,t)dr 



\\ ^\r ,t)^\r\t)V{r-r')^{r\t)^{r ,t)drdr\ 



(1) 



Here V {r-r') is the interaction between the particles, 
which, for the dilute bose gas, is conventionally taken to 
be 



V(r-r')=^oS(r-r'), 



(2) 



Here Uq is related to the scattering length a via 

Uq = ~ — ~ [21]. The potential problem of applying this 

in a self-consistent approach in a 3-dimensional gas is 
described by Huang et al. [22]. In fact our approach 
does not rely on this approximation, as will become 
clear in Sec. 2. For the point of discussion, let us 
however assume, for the moment, that it is valid; this 
will enable us to give our initial discussion in terms that 
link easily with previous accounts. 

The operators ^(r , t) satisfy the equal-time boson 
commutation relations 



[^(r,0, nr\t)] = [^Hr,tl^Hr\t)]=0 



[^(r,0, ^U''',0]=S(r-r'). 



(3) 



(4) 
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The Heisenberg equation of motion for ^(r , t ) then 
becomes 



dt 



2m 



+ V„ap(r) 



+ NUo'F^{r,t)^{r,t) 



nr,t) 



(5) 



This exact operator equation is the starting point for 
most quantum-mechanical treatments of a Bose system. 
We are, in the main, interested in investigating the 
temporal variation of the coherent (mean-field) relative 
to the incoherent parts of the boson field ofi^erator. 
We thus follow the usual route of decomposing ^(r , t ) 
as [6] 



^(r,0=^(r,0 + S(r,0 



(6) 



Here 0(r ,t) = ('^(r, t)) expresses the non- vanishing 
macroscopic value of the Bose field associated with 
broken symmetry due to the presence of BEC. We note 
that there are some subtle extra issues that arise in this 
procedure in the case of small trapped assemblies of 
atoms. For the point of view of most of this article we 
shall assume that this is a well-defined quantity. The 
reader should however bear in mind that in small 
systems <P(r,t) may evolve in amplitude and phase 
[23], [24]. In this paper we shall develop a method for 
determining this evolution. 

Using this decomposition we note that we can re-ex- 
press the interaction term ^^ ^^ of Eq. (5) in the form 



|0|'0-F 2|0|' 8 -F 208^8 -F0' 8^ -F 0*88 -F 8^88 . 

(7) 

The first term in the above expression corresponds to 
the mean field term which, if it were the only term 
retained, would convert the mean value of Eq. (5) into 
the usual Nonlinear Schrodinger Equation (NLSE) or 
Gross-Pitaevskii equation [25] 






+ NUo\^{rj)Y^{rj), 



(8) 



This is the simplest mean field equation describing the 
atomic assembly and is only strictly valid when all the 
atoms are condensed, e.g., at temperatures close to zero. 
The NLSE has already been studied in some detail along 
with the link to the 7=0 elementary excitations 
(phonons) in isotropic [26] - [29] and anisotropic [30], 
[31] traps. 

The next four terms appearing on the right hand side 
of Eq. (7) represent the interaction between condensate 
and excited states. Let us now discuss their meaning. 
Consider substituting those terms into the operator 
equation of motion Eq. (5) and taking the average of that 
expression. The first two terms in the resulting mean- 
field equation correspond to the direct and exchange 
Hartree-Fock (HE) terms. The next two terms represent 
the Bogoliubov correction. In the Popov approximation 
[13], [14] of the HFB equations, the anomalous averages 
of the form (8 (r , t)h{r , t )) are assumed to be negligi- 
ble, so that the second of the Bogoliubov corrections is 
dropped. This approximation can be shown to give 
accurate results for a homogeneous gas with purely 
repulsive interactions. 

Let^us now substitute Eq. (7) in full (i.e., including 
the ( 8^88 ) term) into the operator Eq. (5) and consider 
its mean value. This gives the exact equation of motion 

+ NUo\<^{r,t)\^<^{r,t) 

+ 0*(r,O<8V,OSV,O)} 

+ NUo$Kr ,t)l{r ,t)l{r ,t)) , 

(9) 

Conventionally, the terms (8^88 ) are handled in the 
self-consistent approximation [13] 

i\r ,08(r ,08(r ,0 - 2 $\r ,08(r ,0)8(r ,0 
+ (8(r,08(r,0)8^(r,0. (10) 

This implies that such terms vanish trivially in the finite 
temperature generalization of the NLSE, given by 
Eq. (9). However, the terms (8^88 ) contain the effects 
of more complex correlations between the particles, 
which are neglected in the HFB approximation. We have 
chosen to analyze the effect of these correlations, thus 
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extending the pure mean-field treatments (e.g., Ref. 
[13]). Such correlations of three particle operators are 
transiently produced during collisions and thus generate 
the collisional feed term for the condensate, as discussed 
in [32], [33]. For the dilute gas, these kinetic effects can 
be dealt with, as discussed below and in [20]. We note 
once again that we have chosen to do the calculation in 
this way, because such treatment gives a full equation of 
motion that is valid outside thermal equilibrium and can 
handle essentially arbitrary initial states. The Keldysh 
non-equilibrium approach [16], [17] should clearly 
produce results entirely equivalent to ours in the dilute 
gas limit. One should however bear in mind that the 
results from two such different formalisms may be diffi- 
cult to reconcile at first sight. 

Let us consider the behavior of triple product aver- 
ages over time scales that are long compared to the 
duration of a typical binary atomic collision. A self- 
consistent evolution may then produce non-vanishing 
mean values of the triples, which would mean they have 
to be considered on an equal footing with the other order 
parameters of the system (i.e., BEC or pairing). Such a 
possibility has already been discussed in [15], where a 
cluster expansion approach is used, instead of the usual 
perturbative expansion in the interatomic potential, for 
the description of dilute degenerate gases with short- 
range repulsive interactions. In this paper we derive a 
self-consistent set of equations that allows for the inclu- 
sion of possible non-vanishing triple product mean 
values, which we shall henceforth refer to as "triplets" 
in the spirit of Laloe [15]. Whether such an order 
parameter can be present in a real system and how 
it might affect the behavior of the system will be the 
subject for future study. 

We shall therefore procede in a generalized mean- 
field approach for which triplets are retained as a possi- 
bility. We are particularly interested in laying the ground 
for the study of such self-consistent solutions and their 
role in the "macroscopic" behavior of a gas. When 
looking at the structure of Eq. (9) it is clear that in order 
to obtain a closed set of equations for the various possi- 
ble order parameters, we shall need the evolution of the 
averages of products of two and three fluctuation opera- 
tors 8^^\r ,0- In our self-consistent coupled equations, 
all operators being averaged over have the same time 
label which we shall not, for the sake of convenience, 
give explicitly in the equations. In expressing the equa- 
tions for such averages we have used suffixes 1 and 2 to 
indicate whether the boson field operators are acting on 
spatial coordinates r or r ' respectively. For the binary 
products, corresponding to what we shall term the 
generalized density matrix, we find the following 
equations of motion 



ot 



A^ 



(vi-vo+(v.g^-v.^^) 



<8!82) 



+ NUa { <&I<8I8|) + 2|<&2|' <8I82) + (81818282)} 



-NUo 



2(&2<8!8|82)+(&2 (818282) 



+ NUo {(a):)^<8i82) + 2|a), 1^(8182) 
-(81818,82)} 



+ NUo 



2(&t(8l8,82)+(&i(8l8l82) 



ih— (8182)= 
at 



(v?+vi)+(y£^-y.g^) 



(11) 



(8182) 



-FiVt/o {^2'<8|8i)-f2|^2|' (8182) + (81818282)} 



+ NU^ 



2^2(818182)+ ^2(818282) 



+ 8^2(^1 + (8282)) 

+ iVt/o{^(8l82) + 2|^i|^(8i82) 

+ (81818182)} 

+ iVt/o [2^1(818182)+ ^1(818182)], 



(12) 



here 8^2 denotes a Kronecker delta. 

These equations are again exact, although fairly use- 
less without any approximations. We note that, apart 
from averages of products of three fluctuation operators, 
terms with products of four such operators appear on the 
right-hand side of the above equations. These represent 
more complex correlations between particles under- 
going collisions and we will make appropriate approxi- 
mations for them as we proceed. 

In principle, one could work out the equations for the 
rate of change of both these quantities which would in 
turn depend on higher order correlations. It is, however, 
essential to terminate this coupling to higher and higher 
order parameters by making some suitable truncation, 
which is a very natural procedure in a weakly-coupled 
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gas [7]. In this paper, we have chosen to include the 
possibility of mean values of triple products in our 
model and reduce higher particle correlations to lower 
ones (via an approximate version of Wick's theorem). 
By including all correlations up to the triplets, we can 
describe the effects of the kinetics on the evolution of 
the condensate in a fully self-consistent mean-field ap- 
proach. In order to obtain the desired closed set of 
coupled equations within this framework, we need of 
course to consider the equations of motion for averages 
of products of three fluctuation operators, which is done 
in Sec. 4. 

The method described in this paper can, in principle, 
be used for calculating the relative importance of these 
various order parameters. The possible appearance of a 
pairing or other order parameters in the case of attrac- 
tive interactions has been anticipated for some time [34]. 
In [35] Stoof has discussed the possibility of the anoma- 
lous average (88) becoming the important order 
parameter in the case of a homogeneous gas with attrac- 
tive interactions. We would like to examine the closely 
related issue for a trapped inhomogeneous gas. 

In the next section we repeat the above analysis in the 
occupation number formalism and discuss the various 
approximations made in more detail. We then show how 
the equations can be cast in a generalized density matrix 
form (Sec. 3) and derive equations of motion for the 
triplets (Sec. 4). We further illustrate the consistency of 
the derived equations with simpler ones found in the 
literature (Sec. 5) and argue how our analysis can be 
used to set the basis for a full non-equilibrium theory. 



Here the ipi are the spatial parts of the Fock space 
decomposition of ^(r, t) given by 



hr,t) = ^Mr)Mt). 



(15) 



Analogous to Eq. (6) we now decompose the time- 
dependent single-particle operators according to 



d,(t)= z,(t)-\-c,(t). 



(16) 



This decomposition is extremely useful as it allows the 
application of Wick's theorem for decomposing higher 
order correlations into lower order products. 

Initially we consider the temporal dependence of the 
mean field which yields the exact equation ^ [analogous 
to Eq. (9)] 



2pms{t)Zr{t) -H zl{t)Zm{t)Zr{t) + Kmr{t)zl{t) 



+ {cl{t)c,{t)Crn{t)) 



(17) 



Here we have defined the one-body density matrix 
p{t) and the two-body or pair density matrix K{t) by 
their respective matrix elements^ 



2. Mean Field Equations for a Trapped 
Atomic System 

In the occupation number representation the hamilto- 
nian for pairwise interactions has the form 



//=2 {r\a\s)dl(t)d,(t) 

rs 

+ 2 2 {rs\V\mn)dl(t)dl(t)dn(t)d^(t) . (13) 



In this equation the operator H contains the kinetic 
energy as well as the trap potential and {rs\V\mn) 
represents the interaction potential between a pair of 
particles. This is defined by 



P,{t)^{i\p{t)\j)^{c]{t)c,{t)) (18) 

K,(t)^{ij\K(t))^{c,(t)c,(t)). (19) 

We note that p{t) and K{t) are both nX n matrices, 
where n is the number of single-particle states. 

The equation of motion for the single-particle density 
matrix Py(0 is given by 



i/i^^^=-{[H,a]{t)a,{t)]) 



At 



.r/iz,-^-r/iz,^^ (20) 



(ry |y|mf7) = 



// 



.a; (r).A;(r')y(r -/•')./.„ (r').A„(r)drd/-'. 

(14) 



^ Due to the extensive use of the operators a/^\ c,-^^^ we have dropped 
the explicit operator "hat" notation f). 

^ We would like to point out here that these matrices transform differ- 
ently under a unitary transformation, which is implicit in the notation 
{i\p{t)\j) and {ij\K{t)). As a matter of fact, K{t), unlike p{t), 
transforms as a matrix only under real orthogonal transformations. 
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with 



[H,a](t)ai(t)] = ^ {r\S\j)aHt)a.(t) 

r 

-2 {i\3\r)a]{t)ar{t) 



+ 2^ {rs\V\m])al {t)al {t)a^{t)ai{t) 



Zj {ir\V\mn)a] {t)al{t)an(t)afn(t) , 



(21) 



We now proceed to decompose operators ap(t) accord- 
ing to Eq. (16). Hence (upon suppressing temporal de- 
pendence) the exact form of Eq. (20) contains terms 
with the following structure (with the appropriate in- 
dices): 

(H)[z*z,(c^c)] 

and 

{V)[z*z*zz,{c^c)z''z, {cc)z*z*, {c^cc)z*, (c'^c^cc)] 
and their hermitian conjugate quantities. At this point 
we will simplify our expressions by assuming that quan- 
tities of the form (c ^c ^cc ) can be approximated by their 
mean-field values according to the conventional form of 
Wick's theorem 

(cjc^ ) {c\Cn ) + {c\Cn ) {clc^ ) + {c\c\){fmCn) 

In this approximation, we are essentially assuming that 
the effects of higher order correlations (such as 
(c^c^cW)) on the mean values of (c^cW) can be 
neglected. Within our mean-field approximation defined 
by Eq. (22), we obtain the following equations of motion 



dt 



+ 2 (ir\V\mn)\(c]clc^) Zn 

rmn ^ 

+ {c]clc^) ZmHc]Cr,C^) Z^ 

-2 {rs\V\mj)Uclc^Ci) zl 

rsm ^ 

+ {cUmCi) zl + {clclCi) Z 



Here we have defined the time-dependent Hartree-Fock 
hamiltonian h{t) and the pairing field A{t) by their 
respective matrix elements [11] as follows: 

h.j^{i\h\j)^{i\^\j)+2Z {ik\V\lj)Vzlz&pik^ 

kl 

(24) 
Atj ^ {ij \A) ^ E {ij\V\kl)iKu+ZkZi^ . (25) 

kl 

We note that in an alternative (variational) approach [11] 
(see Appendix A), these may be defined as 



hij (0 ^ 



dE 
dpjiit) 



4,(0 = 



dE 



(26) 



where E = E[z, p, k^ is the energy functional of the 
system. 

The equations of motion for the pair matrix k can be 
derived along the same lines as above. Consideration of 
the term ([//, Uj (t)ai (01) gives rise to the correlation 
{ajalunam) (where r, m, n represent indices that are 
being summed over). By using the decomposition equa- 
tion [Eq. (16)] for obtaining correlations of fluctuation 
operators and imposing the boson commutation rules for 
the c -operators, we obtain the term {c IcjCnCm ), which is 
handled in the mean-field approximation in an 
analogous fashion to Eq. (22). 

Hence, the equation analogous to Eq. (12) reads 



dr 



+ 2 iir\y\mn)\ic\cfn)Zm 

rmn ^ 

-F (flcf^) Zn+{CjCnC^)zl \ 

+ 2 ijs\V\mn)\{clcnCi) Zn, 

smn ^ 

+ {cU^Ci) Zn+{CnCmCi)zl \ ■ 



(27) 



(23) 



The parameters z, p, k and {c^^^cc) can be treated 
(self-consistently) as suitable order parameters for the 
system, provided they are well-defined. This could be 
true either in a completely static or a slowly-varying 
situation. It is however also possible that the kinetic 
correlation terms {c^^^cc) are established on time-scales 
much shorter than the mean evolution times of the 
parameters z,p, and k . This is equivalent to saying that 
the collision duration relevant to the evolution of 
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{c^^^cc) must be much smaller than the mean time and the matrix i] by 

between collisions. This should hold for dilute weakly- 
interacting Bose gases for which n^d? « 1 where n^ is > 
the number density in the trap. _ / 
In the next section we write these equations in a more \ 
compact generalized matrix form. We then derive equa- 
tions of motion for the triplets (Sec. 4) and show that if 
the latter are neglected (Sec. 5), the equations derived in 
this paper are consistent with simplified versions 
appearing in the literature. 



(34) 



3. The Equations of Motion in a 
Generalized Density Matrix Form 

Equations (23) and (27) can be cast in a generalized 
matrix form 



i^^= [/?,p]-(/czl*-zl/c*)+(M-M*) (28) 



where U is the unity matrix of the fz -dimensional Fock 
space. 

For convenience, we have also defined the boson 
quasiparticle hamiltonian 



H-- 



(35) 



We note that the matrices R and H are not the same as 
those defined by Blaizot and Ripka [11] as discussed in 
Appendix A. 

We have further defined the matrix 



i^^= {hK+Kh'')+{pA+Ap')+A+{N+N) (29) 



where the matrices M and A^ have been defined in terms 
of their elements 

^y = 2 {is\y\mr)i{c]clc^) Zr 

rms 
+ (cJcUr) ZmHc]CrC^) ZU (30) 

^y = 2 {is\V\mr)[{cUmCj)zr 

rms 
+ (ctCrCj) Zm-^{CrCmCj) ZU (31) 

and N represents the transpose of matrix A^. 

It is straightforward to show that these can be written 
in matrix notation as 



d/? 



(32) 



where all matrices are of the form 2nX2n. Here we 
have defined the generalized boson density matrix R by 



K-- 



M-M" (7V+7V) 



(36) 



which gives the effect of the triplets on the single- 
particle and pair excitations. 

Thus far, we have derived equations of motions for p 
and K which involve triplets. In order to get a closed set 
of mean values, we need the equations of motion for the 
triplets which we derive below. 



4. Equations of Motion for the Triplets 

Before deriving the equations of motion for triplets of 
the form (c^'^^cc) we would like to point out that they 
can be used in two ways. Firstly, they allow us to con- 
struct a self-consistent generalized mean-field theory. 
Furthermore, they offer the possibility of including 
kinetic effects of condensate formation. 

We define the time-dependent nXnXn rank 3 triplet 
tensors y(t) and \(t) by their respective elements 



yijk = {ciCjCk) 



(37) 



and 



R = 



-k' Ap'+a)^ 



(33) 



Xijk = {c]cjCk) 



(38) 
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The equations of motion for the triplets can be ob- 
tained by following the procedure of Sec. 2. We note that 
this leads to very complex expressions, so the account 
will be kept very brief here. 

Consider initially the equation of motion fot the 
triplet jijk . Upon imposing the decomposition equation 
[Eq. (16)], we note that the right hand side of the result- 
ing equation contains — apart from the simpler correla- 
tions — correlations of four and five single-particle fluc- 
tuation operators. We assume that the effects of 
correlations of six or more operators (present on the 
right hand side of higher equations of motion) are negli- 
gible when looking at the behavior of the system during 
a typical collision. This means that we set all terms 
(cccccc) — ({cc)y, {{cccyf (appearing in the equations 
of motion) to zero, irrespective of their indices or 
whether they contain products of normal or anomalous 
averages. Correlations of four fluctuation operators are 
handled in the mean-field approximation Eq. (22). In 
treating the equation of motion for y we also use Wick's 
theorem to decompose correlations of five fluctuation 
operators in the form 



\C s C iCj^CfiCffi/ Pis /knm ' Pksjinm "■" Pns jikm "■" Pms jikn 
+ KkiAsnm "T f^ni^skm "•" ^mi^skn "•" ^nk^sim 
+ f<mk^sin + f<mn^sik ■ (39) 



It is essential that the above decomposition is in accord 
with the mean field approximation defined by Eq. (22). 
Strictly speaking Eq. (22) should contain an extra term 
on the right hand side due to the effects of five particle 
operator correlations, and these would have to be in- 
cluded in a treatment that went beyond the mean-field 
approach. 

Writing the resulting equation of motion in more 
compact notation we obtain 



/^ ^ ( TyO = [h,y},jk-\-{A,\} ijk-^(W,jk-\- Wjk,-\- Wkij ) 



+(A.+r,.,+A,) (40) 

where we have defined the following quantities 

[h, y}ijk= {hy)ijk+{hy)jki+{hy)kij 

= 2 (^ir yrjk + hjr yrki + kf^r 7ry ) (41) 



= 2 i^irXrjk + AjrXrki + AkrXrij) (42) 

r 

^ijk = 2 {is\V\mr)[KnijLskr-^ KmkLsjr 
rms 

+ Pjsykrm+ Pksyjrm] (43) 



^rms — ^PrmZs "•" Zr ^ms "•" ^^r. 



(44) 



^ijk = 2 W \yW^)l'^f<kmZs-^ykms] ■ (45) 



Similarly we have 



i^— (Xijk) = (Xh*)ijk-\-(Xh'')ikj-(h*X)ijk 



+ (A*zl),,,+ (A*zl),,^-(zl» 



ijk 



-F (Xkj,-\-Xjki - Y^k )-^Ajki (46) 

^kji = 2 (^^ \y\mr)[K^jJris + PriLsjm 

rms 

+ K-syrmj-^ PjsXmr] (47) 

Yijk = 2 (^^ \y\mr)[p^jJrks + PmkJrjs 
rms 

-^ K*jXk^,-\- K^kXj^,] (48) 

J rms — -^Zr^ms "r PrmZs "r -^^rms \^^ ) 

Ajki = 2 (jk\V\ms)[2p^iZs-^Xi^s] . (50) 
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There are a range of ways in which we can make use 
of these equations. Firstly, we could solve Eqs. (17), 
(23), (27), (40), and (46) self-consistently as a general- 
ization of the HFB equations that includes the extra 
anomalous averages. Exploring the phase diagram for 
the various order parameters will, however, in itself be 
an enormous task. It will, nonetheless, be an essential 
one in the study of these new evaporatively cooled 
assemblies. In the first instance an examination of the 
HFB-Popov behaviour will be an important advance. 
The next step should be the study of how {cc) and 
(c^^^cc) arise in the presence of attractive interactions. 
We believe that this will be an important theme in the 
next few years. We can also use these equations to inves- 
tigate the damping of time dependent solutions that 
come from one-loop processes, e.g., the Landau damp- 
ing discussed by Payne and Griffin [36] for the homoge- 
nous case. This means we can explore the range of 
validity of the simpler NLSE approaches to condensate 
time evolution. 

These equations can also be used to treat condensa- 
tion kinetics. It should however be pointed out that a 
completely general treatment that also includes kinetics 
of excited states, further requires us to look at the equa- 
tions of motion of terms {c ^^^c ^^Vc ), where the 
algebraic complexity becomes quite formidable. 



5. Reduction to Simpler Equations 

We now neglect the triplets and obtain simpler equa- 
tions in appropriate limits, some of which can be found 
in the literature. 

5.1 The Nonlinear Schrodinger Equation 

If we further ignore single-particle p and pair k exci- 
tations in the time-dependent equation for Zi{t) [Eq. 
(17)] and restore explicit time-dependence we end up 
with the mean-field equation 



ifi^^=J.hfmZr{t) 



(51) 



where the time-dependent mean-field Hartree-Fock 
hamiltonian h ^^\t) is defined by 



h^Sit) ={i\S\r)-^{is\V\mr)z:(t)Zm(t) . 



(52) 



One should note that the interaction potential appearing 
in this equation is the original one. The replacement of 
this by the two-body T-matrix is straightforward but 
relies on a decoupling — ladder approximation — in the 
equations of motion. We shall give the explicit account 



of this in a future publication. With this replacement 
Eq. (51) becomes the analogue of the NLSE [Eq. (8)]. 

0(r, 0+ A^t/o|^(r, Ol' ^(r, (53) 

in the occupation number representation we note that 
Eq. (15) can be equivalently written in the form 



0(r,f) = 2<^^(r)zKO 



(54) 



This confirms that we can reconstruct the NLSE by 



multiplying Eq. (51) by 



Vn 



and summing over all 



indices /. The factor Vn has been included so as to 
re-scale the time-dependent coefficients Zi(t) which 
obey the normalization condition 



^zkt)Zi(t) = N . 



(55) 



At r= the whole system is condensed in the lowest 
energy eigenstate and since according to Eq. (54) all the 
time-dependence has been included in the coefficients 
Zi (t), we can re-express these as 



i/it 
Zi(t) = ZiC^ 



(56) 



where the Zi are now simply constant complex numbers 
and fjL is the ground state chemical potential. Clearly 
Eq. (51) then reduces to the time-independent NLSE 

[11] 



^hfhr= fJ^Zi 



(57) 



Here we see that we can trivially transform from the 
time-dependent to the time-independent picture, merely 
by setting the time-dependent term on the left-hand side 
to zero and simultaneously replacing the time-depen- 
dent mean-field Hartree-Fock hamiltonian 

hf(t)^{i\S\j)^^ {ik\V\lj)zl(t)zdt) (58) 

ki 

by the time-independent one given by 

hf^{i\S-fi\j)^^ {ik\V\lj)zlzi . (59) 

ki 

We note that this remains true even at non-zero temper- 
atures, when the second term also contains correlations 
due to the excited particles. 
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5.2 Static Equations 

In this section we will show that Eqs. (23) and (27) 
are also consistent with the r> static HFB equations 
that we have additionally derived from a variational 
principle approach in Appendix A. Neglecting all time 
dependence in Eqs. (28) and (29), we thus obtain the 
static matrix equations 

[h,p]={KA''-AK) (60) 

[// K + Kh"] + (pA + Ap*) + Zi = (61) 

However, to solve these equations self-consistently, 
we also need an equation normalizing the vector z. 

We use the substitution Zi{t) = Zi^ ^ in Eq. (17) to 
obtain the static equation 

^{i\E-p.\r)zr+ 

r 

2 {is\V\mr){2p^sZr+zlzmZr-^K^rzl} =0 (62) 

rms 

where the chemical potential [x appears as discussed 
earlier. 

This can be re-expressed as 



^{KZr^AfrZ:) = 



(63) 



where 



/i5^(/|H-)a|7) + S {ik\y\lj)Vzlzi+ 2pi,] (64) 

ki 

and 

A[j^^{ij\V\lkl)Ku. (65) 

ki 

More generally, we obtain the matrix equations [11] 

(66) 




which give us the values of all Zi for fixed values of h ^ 
and A"". Equations (60), (61) and (66) constitute the full 
set of static Hartree-Fock-Bogoliubov equations which 
can be solved self-consistently. If we set Zi = /c = 0, we 
thus obtain the static self-consistent HFB description of 
a condensate with elementary excitations in the Popov 
approximation [14], which further reduces to the time- 
independent NLSE for p = 0. 



6. Discussion 

In this paper we have developed a time-dependent 
self-consistent generalized mean-field method for 
describing a trapped and partially condensed atomic 
system. Central to our approach lies the operator decom- 
position equation [Eq. (6)] (or equivalently Eq. (16)) into 
a mean field and a fluctuating part. This decomposition 
is reasonable, as we do not make any assumptions about 
the initial value of the wavefunction (r ,t), or its 
evolution. Consequently, our treatment should remain 
valid over a wide range of conditions, e.g., for cases 
close to the transition temperature Tc. Our equations 
allow for the possibility of other order parameters and 
give us a way of investigating their possible role in an 
evaporatively cooled gas. 

In our treatment, we have considered an inhomo- 
geneous gas, which implies that the excited states of the 
system are strongly influenced by the finite size of the 
trapped assembly. This natural length scale will, of 
course, lead to the presence of a gap in the excitation 
spectrum, as opposed to the gapless spectrum of the 
homogeneous gas. 

Let us initially assume that the triplets are small in 
comparison with the usual HFB order parameters. In the 
case of moderately- sized condensates, the solutions of 
the HFB equations will hence be a good basis, as long 
as the kinetic processes produce collisional widths that 
are less than the separation of the levels. However, as the 
number of atoms in the trap is increased, we will 
approach the so-called quasi-homogeneous limit for the 
condensate. This happens when the healing length for 
the gas is much smaller than the size of the condensate. 
In that regime, it becomes essential to describe the mean 
fields by means of local densities [i.e., z(r), p(r), and 
/c(r)]. We thus obtain locally well-defined quasiparti- 
cles (in the sense that their collisional widths are smaller 
than the characteristic frequencies) in regions specified 
by the healing length of the atomic assembly under 
consideration. This approach has been employed by 
Dorfman and Kirkpatrick to derive transport coeffi- 
cients for the dilute Bose gas [32]. 

The discussion above relies on assumption that the 
HFB states for the trap as a whole constitute the most 
sensible state. This, in turn, implies that we should, in 
order to be consistent, assume that any kinetic rates are 
less than the frequency separation of the HFB excita- 
tions, which seems reasonable for the case of conden- 
sates with modest numbers of particles. Our equations 
are, however, not limited to the pure HFB regime. For 
example, in the case of repulsive interactions, both pair- 
ing and triplets may be treated purturbatively and the 
analysis extended beyond HFB. 
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7. Appendix A. Proof of the Static 
Generalized Equations using a 
Variational Principle 

In this Appendix we derive the static form of the 
generalized quasiparticle equations of motion using the 
variational principle approach employed, in detail, by 
Blaizot and Ripka for fermions [11]. 

Consider a system in thermal equilibrium described 
by a statistical density matrix D given by 



which can be summarized into 



Tr(DQ)=-tr(HR) 



(74) 



where we discriminate between the trace operation in 
Fock space (7r ) and in the space of single-particle states 
(tr). Thus, 



8[ln(Zo)]=-§rr(/?8H) 



(75) 



D- 



Zo 



(67) 



where Zo is given by 



Zo=7>(e-'*«) 



(68) 



Here R represents the generalized density matrix of a 
quasiparticle vacuum and H the quasiparticle hamilto- 
nian, respectively defined in Eqs. (33) and (35). 

We note that these are not the same as the ones 
defined in [11], but they are related by 



R = vRo 



H = Hor] 



(76) 



and Q is the quasiparticle operator 

Q=o^ {hij(cJcj-^CjcJ)-^(AijcJc]-^AijCjCi)} (69) 

whose matrix elements hij and Aij (defined as in Eqs. 
(24)-(25)) will be considered as variational parameters. 
In applying the variational principle, we seek to mini- 
mize the thermodynamic potential 



where the zero suffix has been used to represent the 
quantities defined in [11]. The matrix 7] appears here 
due to the nature of the boson commutation relations 
that must be satisfied by the transformed quasiparticle 
operators. 

We also note that, for a particular state (which can be 
explicitly constructed from the particle vacuum), the 
generalized density matrix R obtains the canonical form 



(I)(D) = E-ijlN-TS 



(70) 



where E=Tr (DH) and N=Tr (DN) are respectively the 
system's average energy and particle number (H has 
been defined in Eq. (13) and N = XiaJai is the total 
number operator), Tis the temperature and S the entropy 
of the system. 

The functional we shall minimize is thus 



cj> = -^\n(Zo)-Tr(DQ)^ E 



(71) 



R = 




^0 -7?. 



(77) 



which satisfies the conditions 

R'^ = R R^ = -R 



(78) 



Using Wick's theorem for ensemble averages, we can 
also express the energy E in terms of /?, so that the 
minimization of the functional cj) 



which represents an upper limit to the actual partition 
function Z. Here E = Tr {DH ^) is the expectation value 
of the hamiltonian H^ = H -fiN . 

We note that the single-particle (p) and pair density 
matrices (k) are given in the statistical density matrix 
representation of a system by 



gives 



8(A = irr(//8/?) + ||a/? = 



V-¥-H 



(79) 



(80) 



p,j = Tr[Dic]c,)] 



(72) where we have used the notation [11] 



Kij=Tr[D(cjCi)] 



(73) 



tr 



1^)H-| 






(81) 
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We hence deduce that cf) is stationary with respect to 
arbitrary variations of H (or equivalently of R related to 
H)if 



dE[R] 1 



(82) 



At steady state the energy functional E [/?] remains 
constant, so that the second constraint of Eq. (78) can be 
taken into account by introducing a matrix of Lagrange 
multipliers L as 



8 {£[/?]-rr[L(/?' + /?)]}=0, 



(83) 



Eliminating the Lagrange multipliers, we thus obtain 
the equation 



[//,/?] = 



(84) 



which simply expresses Eqs. (60) and (61) in a more 
compact form. In order to obtain the full set of the 
steady state Hartree-Fock-Bogoliubov equations, we 
also need an equation for the vectors z. This is obtained 
in [11] by consideration of the variation of the energy 
functional E \R\ with respect to z and z*. The resulting 
equation is identical to Eq. {^6). 
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